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FOREWORD 


This document is the final report of work done on adaptive finite 
element methods for SSME internal flow analysis for NASA under Contract 
NAS8-36647. The results described here summarize a one-year pilot study 
into several classes of adaptive methods which may have important implica- 
tions in a large body of computational fluid dynamics work connected with 
the analysis and design of the space shuttle main engine. The one-year 
effort supported by this NASA-Marshall Space Flight Center contract repre- 
sents one component of a larger program on adaptive methods in CFD underway 
at the Computational Mechanics Company, Inc. of Austin, Texas. The very 
encouraging results obtained suggest that a significant increase in the 
reliability of computer-generated flow simulations is possible through 
the use of special error estimates and various adaptive procedures. 



1. INTRODUCTION 


In this report, we describe adaptive finite element methods for the 
analysis of classes of problems in compressible and incompressible flow 
of interest in SSME (space shuttle main engine) analysis and design. 

The general objective of the adaptive methods of interest here is to 
improve and to quantify the quality of numerical solutions to the governing 
partial differential equations of fluid dynamics in two-dimensional cases. 
Implicit in this goal is the resolution of two questions: 1) how is this 
quality of the numerical solution to be assessed? and 2) how does one 
effectively adapt the solution to improve quality? 

An answer to the first question is to somehow estimate the local 
error in a given mesh cell (in an appropriate norm) by using the results 
of a first trial calculation. Stated in another way, one can construct 
a-posteriori error estimates using data from rougher solutions obtained 
on a coarse mesh or with a low-order approximation. Once an estimate of 
the local quality of the solution is known, one can "adapt" the numerical 
scheme so as to improve the quality of the solution; i.e. so as to reduce 
the local approximation error. 

There are several different families of adaptive schemes that can be 
used to improve the quality of solutions in complex flow simulations. 
Among these are 1) r-methods (node-redistribution or moving mesh methods) 
in which a fixed number of nodal points is allowed to migrate to points 
in the mesh where high error is detected; 2) h-methods, in which the mesh 
size h is automatically refined to reduce local error and 3) p-methods, 
in which the local degree p of the finite element approximation is increased 
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to reduce local error. There are, of course, combinations of these methods 
which can be very effective. Two of the three basic techniques have been 
studied in the project reported here: an r-method for steady Euler equations 
in two dimensions and a p-method for transient, laminar, viscous incom- 
pressible flow. For discussions of our work on h-methods for these types 
of problems see, e.g. [ 1, 2, 3, 4]. 

The issue of a-posteriori error estimation is a difficult one. Two 
basic categories of error estimation were studied in the work described 
here: 1) residual methods and 2) interpolation methods. The former classes 
of methods make use of residuals computed using trial finite element 
solutions. These methods can be computationally expensive; however, they 
can yield very good estimates of the local error. Interpolation methods, 
on the other hand, are easily implemented but may yield quite crude 
estimates of the actual error. These schemes employ interpolation theory 
and exploit superconvergence properties of finite element methods; they 
seem to be perfectly adequate as a basis for adaptive mesh schemes for 
the classes of problems considered here. 

Following this Introduction, we present weak forms of Navier-Stokes 
equations and the Euler equations that are used as a basis for the 
development of finite element approximations. Since we anticipate the use 
of techniques which move nodes and elements, we construct general space-time 
"variational formulations" for these problems for which the computational 
domain can vary with time. In Section 3, we present two types of adaptive 
schemes. First, an r-method for two-dimensional steady problems in inviscid 
compressible flow characterized by the Euler equations. Numerical results 
are presented for some representative test problems. In Section 3, we 
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also outline a p-method for incompressible viscous flows and cite some 
preliminary numerical results. A brief introduction to residual methods 
of a-posteriori error estimation is given in an Appendix, and some pertinent 
conclusions of the study are listed in Section 4. 
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2. SPACE-TIME VARIATIONAL FORMULATIONS 


OF COMPLEX FLOW PROBLEMS 


2.1. A Space-Time Navier-Stokes Formulation . A general space-time 
variational principle for incompressible viscous flow is characterized as 
follows (see [ 1, 2, 3]). 

Find a velocity field u in a class of functions V such that over a 
time interval [ 0, T] , 

{ p [ v t , u E ] t + p(( u E , v )) t + pb ( u E , u E , v ) 

0 

+ e 1 ( div u E , div v ) } dt 

T 

( f, v) t dt + p (u, v) Q - p (u E , v ) T 

V v e V (2.1) 




where v is an arbitrary test function, p the mass density, M the 
viscosity, f the body force, and 


[ 


V 



(( u, v)) t = 


( u, v) 



u • v dx 


u dx at time t 


V v dx at time t 


at time t 


n 


t 


b (u, 
b (u, 


u, v) = b (u, u. 


v, w) = [(u *V ) 


v) + b (u, ip, v) + b (<J>, u, v) 
v • w + k. div u (v • w) ] dx 


'fi 


t 
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with the spatial time domain at time t. Here we use a penalty method 
( artificial compressibility ) to approximate the hydrostatic pressure 
P, by 


p = - e div u 
e e 

with e a small positive parameter. The functional b (•,*,•) represents 
the convective term in the Navier-Stokes equations and ij) is a particular 
function designed to simplify the enforcement of no-flow boundary con- 
ditions . 

A finite element approximation (2.1) is obtained by replacing 
u E and v with appropriate discrete approximations defined over a space- 
time element K; e.g. 


^ ® “h (x ’ (X) 


N 

Once a finite element solution is obtained on a fixed mesh, we use 

it to compute a local error indicator <p„ which bounds the local e* 1 = u E - 

” K 

in an appropriate norm: 


K 


I 4 


K 


for element K 


We shall discuss means for obtaining <j> 


K 


later. 


2.2 A Space-Time Variational Formulation of the Euler Equations . 
By following a plan similar to that used in the formulation of the 
weak-space-time problem (2.1), a space-time formulation of the Euler 
equations in two dimensions can be obtained. 
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If U(x,t), (x,t) C D, is the 4-vector of conservation variables, 

T 

U = {p, m, E} , with p the mass density, m the linear momentum, and 
E the total energy, and if d£2 and dS denote Lebesque measures of area 
(volume) and length (area) of £2 and 3£2 respectively, then we demand 
that U satisfy the following system of conservation laws: 

UdS=- Q(U)n dS (3.1) 

n~ ha 


d_ 

dt 


Here, Q(U) is the flux and n is the unit outward normal to 3fi . If 
m^ , denote Cartesian components of m , then 

T 

U = {p , mj , m 2 , E} 


Q(U) = 


m l 

' m 2 

p _1 m^ + p(U) 

1 - 1 
' p ni-m, 

i 1 

-1 

-i 2 . 

p m^m 2 

| P m 2 + 1 

p“ 1 m 1 (E + p(U)) 

| p" lm 2 ( E + 1 


n = (n 1 , n 2 > T ; p(U) = (y - 1) (E - p 1 m • m/2) 


In these equations, p is the thermodynamic pressure and y is the ratio 
of specific heats, assumed here to be constant. In addition to (3.1), 
D must satisfy an entropy production inequality as well as an initial 
condition, 

u(x, o) = y Q (x) , x g n 

where Uq is given. 
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It is of fundamental importance to note the smoothness requirements on 
U in order that (3.1) make sense mathematically. Conservation laws (3.1) 
hold when the components of U are bounded measurable (with respect to 
Lebesque measure in x ) functions on D . Thus, we may seek solutions in 
the function space 

V = {V = {V 1§ v 2 , v 3 , v 4 ) t I v ± = v i (x,t) 

G L°°(0,T ; L 1 ^)} ; i = 1, 2, 3, 4} 

In particular, (3.1) is not equivalent to the classical Euler equations, 

H t + div Q(U) = 0 (with U = 3U/3t and div Q = | 3Q al /3x i^ slnce 

solutions may not possess derivatives across surfaces in D . However, 
the conservation laws and initial conditions are fully equivalent to 
the following weak boundary-initial value problem: 


Find U 6 V such that 


m 

(u <L + g(u) : v^)dndt 

D 


V T J(-,0)dfl = f i F T A dS dt 


for all £ G W 

where F is the actual prescribed flux through 3fl and W is a suitable 
space of test functions- 


Here, we use the notation 
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Consider an arbitrary time interval [t^»T 21 d [0,T] and include 
in W functions <j>(x,T 2 ) x 0 • bet u be a subset of fi such that 
(jo n vj T, = 0 > and let F = Qn . Then another weak statement of the system 

k k ^ 

conservation laws over tu x [t^»T 21 is: 

Find U 6 V t0,T such that 
t 2 

(-u T it + (div ^ dn dt 

t j ■'a) 

+ (U T (.,t 2 ) - H T( * » t i) $(*»Tj))dfl = 0 (2.2) 

Ju> 

for all $ 6 

with v“’ T and V U,T appropriate spaces of trial and test functions. 
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3. ADAPTIVE SCHEMES 


3.1 Finite Element Approximations . A hierarchal finite element 

method designed for use in a p-method is used to construct approximations 

of (2.1). We furnish some details late in this section. For a more 

complete discussion of our approach, see [5]. 

Turning to (2.2), finite element approximations of the gas dynamics 

problem are obtained by a direct approximation on finite-dimensional spaces 

approximating the spaces V and W . The spatial domain D is partitioned 

into a collection T, of finite elements H over which the components of 

h e 

trial functions V are approximated by polynomials of degree k . In this 
way, we construct a family (V, ) of finite dimensional spaces of the type 

\ = {y h = {V*. V$, V*. vJ} T 6 V I 6 T k oy , i - 1, 2, 3, 4} 

where P^Oy is the space of polynomials of degree k defined over . 

Alternatively, we can use € Q k ( n g ) » where Q k (y is the space of 

e 

tensor products of polynomials of degree k on ft fe.g., Q, (0 ) is 

e v 1 e 

spanned by bilinear functions, (^(O ) *>y biquadratics, etc.) . In addi- 

tion, a family (W^) of finite dimensional spaces of test functions is 
also constructed. We then consider Galerkin approximations by seeking 
solutions in v h , with W replaced by wj 1 . 

We next derive a special semi-discrete, weak formulation which 
provides the basis for the construction of a popular family of finite 
element schemes. We proceed with the following steps: 
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i) Partition the time interval [0,T] according to 0 = t^ < t^ 
< t 2<***< t N = T; 

ii) Apply the weak balance law (2. 2 ) to a typical time interval 

[ V W ( ” lth T 1 * C n and t 2 ■ C n+1 ); 

iii) Set <j> = 0 in (2. 2 ) suggesting the ultimate use of a 

time-invariant grid (we relax this assumption later); 

iv) Replace the time integrations in (2.2) by the elementary 
midpoint quadrature rule 


n+1 


f(t)dt ~ At f 


n+4 


n 


At = t. - t , f - = f(t + At/2) 
n+I n n 


Thus, with a) = S2 , we obtain the semidiscrete approximation 


& r l * - 


u n dn + At 
n J n 


I 9 

j o 


n+*5 

I • 


: v$ h dn 


- At I tltf* s) d S 

hn 


for all ^ 


(3.1) 


where = U^(x, t^) , etc., being the approximation of 
U , and g n+ 

r h - s(Hr % ) 


n' ' ' ~h 

is the flux at the half step. 


v) To obtain an approximation ' , we use (2.1 ) again for time 

interval [t , t , ] , this time replacing the time integrals bv 
n n+*5 

a simple strip rule and integrating by parts the divergence terms 
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(3.2) 


il nr" - f 

4b 

e e 

$J(div g n )dn 
f or all 

We thus arrive at the algorithm, 

1) With (uJJ , g n = 0(u”)) known at the n th time step, 

compute using (3.4) 

2) Compute Q n+l5 using (3.3) 

3) Compute uj| + 2 using (3.4) 

4) Go to 1) 


A_t 

2 


This algorithm is the finite-element based two-step Lax-Wendrof f / 
Taylor Galerkin scheme. It is one of a family of methods advanced by 
Donea [6], studied by Baker and Kim [7], and successfully refined and 
used by Lohner et al. [8,9] in finite-element applications in fluid 
dynamics. This semi-explicit method is of second order in time and 
can experience spurious oscillations near shocks and other types of 
irregularities in the solution. These deficiencies must be reckoned 
with in implementing the method, and for this purpose we append to 
the right-hand side of (3.1) an artificial term of the form 


- At 


r 4 2 

/II c. (u) U n+1 . 

4 a=1 i=i 1 a > 1 0)1 


d SI 


i 


-At 'T <jv (c (u n )«V U n ) n ds 

an 


with 


3u . 

c i< u > = c a^T 


(no sum on i) 
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In this work, we use meshes of four-node quadrilaterals over which 


the components of U are piecewise bilinear functions. Similar approximations 
and algorithms are used by Bey et al. In addition, so-called group 
approximations of the flux Q (a = 1,2, 3, 4; i = 1,2) are employed 
so that these components are also piecewise bilinear functions determined 
by their values at element nodes. In general, this finite element approximation 
will be of the form, 

N 

U a = l U a (t) VJ> 
j-1 3 


•ai 


N 

- I Q j 
j-l 


ai (t) V 2) 


where N denotes the total number of nodes in the discretization, and' 
i i li li 

a ’ ^ai are va l ues °f U Q at node j , and are the global 

piecewise bilinear basis functions. 

As noted earlier, we advance the solution in time in two steps. It is 
important to note that the first step is essentially local, computed over 
each element, while the second is global and contains the artificial 
viscosity terms: 

First Step : For each element , calculate a constant element vector 


U n+! 5 £ rom 

a,e 




U a + S f dn = l f([ M fl ) u a’" 

a ’ e J n i=i n 1 a 

e e 

-r i 4 dn) <i n \ 
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by solving the 


Second Step ! For each node j , calculate 
following system of equations 



+ 


3d>^ 



dfi} u J’ n+1 


N ( 

l ( 


j = l 


4> .4> . dfi) U j,n 
i] ’ a 


■■ I, *• L ■ 


- At 


an 


n 6 Q «B *1 ds 


Here, Q n denotes the elementwise averaged value of the flux. The 
coefficients t d are defined to be constant over each element, 

p 


T 


8 


cA 

e 


au 

ax 


h 

e 

8 


where c is a global constant (c = 1 in the examples) , denotes the 

area of , ujj denote the components of the fluid velocity. 


3.2 A Node Redistribution Method . We now describe one of the principal 
aspects of this investigation: a moving-mesh, node-redistribution method 
based on equidistribution of error indicators. We begin by presenting 
a simpler mathematical justification of the concept of equidistribution 
of error. 

Consider a regular mesh of quadrilateral elements ^ with a diameter 
h g Let 0 g be an error indicator for element e and suppose that the mesh 
contains a fixed number M of elements. Let h = h(x^, x^) be a mesh 
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function such that 


h(x ,x ) = h = dia(fi ) for (x ,x ) G 0. 

1 L e e t i e 

and note that, approximately. 


M = 





with dfi = dx^dx^ (this being exact for domains which are unions of square 

elements). Let 0 = 0(x^,x 2 ) be mesh function which gives the local error 

indicator when evaluated at a point (0=0 for x G Q ) . We wish to 

e ~ e 

minimize the total error indicator functional, 

M r „ 

J(0) = I 0^ dfl 
e=l J n 

subject to the constraint (3.3). Using Lagrange multipliers, this leads to 
the optimality condition. 


6 (j + A (j 



dfi - M)) 


0 , 


or 


I * L 

e J n 


<5hdfi = 0 


or 


- 30 

n 0 - A = 0 

e e 3h 


2 

Suppose that meas(Q ) = a h and that 0 is of the form 

e o e e 

0 = h G f(u) . Then, integrating this last result over a typical element 

e e 

gives 

ah 3 0. h a_1 f (u)dQ = Ao h 2 
J ^ eee oe 

e 


Hence, the optimal mesh size distribution results when 
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Figure 1 . Calculation of area center-of-error 
N 

X to equidistribute element error 
indicators in a cluster of four elements. 
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0 Z dfl = \a h = CONST, 
e o 


(3.4) 


In other words, to obtain the optimal mesh, we must equidistribute the 


indicators 0 

e 


To use this result to redistribute nodes, we proceed as follows (cf. 


Diaz et al. : 

1) Generate an initial (generally regular) mesh with a fixed number 
M of elements and compute a trial solution on this mesh at one time step; 

2) Compute the corresponding error indicators 0 g ; 

3) For a group k of P elements (with P always 4 in this 

work), let A denote the area of element i in the group. The 
e i 

area-weighted indicators for group k are the P-numbers, 

0 ./A . 
ex ei 

4) Let y denote a vector from the origin of a global coordinate 

e i 

system to the centroid of element e^ of group k . Then the center of 
error of group k is defined as the vector (see Fig. 1) 


4 v> 

e ' 

I i e r 

i=l e i ' e^' 


l r 1 

i=l ( 


5) Relocate the node at the center of group k to lie at the vertex 

f k 
of x ; 

6) Continue this sequence of operations over each group h of four 
elements until the new location of each node does not change more than a 
preassigned tolerance. 
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There remains only the issue of how the error indicator 0 can 

e 

be calculated in an efficient manner. Instead of using residual methods 
such as those discussed in the Appendix, we shall use an interpolation 
method. These methods are derived from the theory of interpolation 
of finite elements (see Oden and Carey [10]. In particular, let u be 
a smooth function defined over a regular domain R. The W r,P (^-semi- 
norm of u is defined by 


W 


r,P 


(A) 


1 


I 

A i+j-r 
i.j*0 




u 


Sx l S *2 


1/P 


d0 


where l ^ p ^ °° and 
The Sobolev norm of u 

M r n 

W r,P (fi)u 


r is a non-negative integer, 
is 



Let G be an arbitrary convex subdomain (a finite element) of 0 
over which u is interpolated by a function which contains complete 

piecewise polynomials of degree k . Then, it can be shown that 
the local interpolation error in the W n,p (G) -semi-norm is 


where 



W® ,P (G) 


< 




W^’P 


CG) 


h = the diameter of the domain G 

P = the diameter of the largest sphere that can be inscribed inside G 
n = the dimension of the domain R 
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p’ = p/(p - 1) 

C = a constant independent of h , P , and u . 

If p is proportional to h and if it remains proportional in refinements 
of G defined by parametrically reducing h , we have 


— — + k+I - m 


m,p,G 


SChP P 


u 


k+1 ,p 


(4.13) 


with l*L „ o “ M m „ » etc. and E = u - u, . 

m,p t G V^’^G) n 

Such estimates can be used to devise crude adaptive schemes. Suppose 
that u on the right side of (4.13) is replaced by a finite element 
approximation u^ and that p ~ ^k+1 p + • Then, (3.4) 

indicates that the local error in the W^’^G) seminorm is proportional to 

Some choices are: 


, . .. ,n/p' - n/p + k+1 - m i 

the error indicator, h * y u 


k+1 »P * 


i) n = 2, m = 0, k = 1, p = p ' = 2 


! E II L 2(g) = ^ U |U| 2 2 


c h 2 |u| 2>2 ^ G C0 G 


(3.5) 


In this case, one must approximate the W 2 * 2 -semi-norm of u over G ; 
i.e., the L 2 -norm of second partial derivatives of u . 
ii) n = 2, p = ”, p' = 1, k = 0, m = 0 . 


|E h | i r . = Ch 2 | E h | 

L 1 (G) average 1 


£ Ch 3 |u| , 

1 


= Ch 3 max|v*u(x)| “ 
x6G 


(3.6) 
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3.3 Numerical Experiments - An r-method. We shall now cite some 


representative examples in which the r-method described above is used. 

All examples here are steady-state problems and the following conventions 
are used: 

1. The numerical solution is computed on a fixed mesh and is advanced 
in time until a steady state is reached. 

2. After convergence to a steady state, initial error indicators 
0^ are computed according to 

0 = A f Vp h • Vp h dfi 

6 6 h 

e 

in analogy with (3.5). 

3. Then, a modified error indicator 0 g is employed which is designed 
to be always greater than unity even when 0^= 0. In particular, we 

use 

a9 e 

®e " 1 + 6 + Y© e 

In our examples, a = 81, 3 = 1, and y = 8. 

4. Nodes are redistributed a total of K times using the procedure 
described earlier. In the example, we take only two iterations (K = 2). 

A Shock Reflection Problem . We begin with a problem for which 
an exact solution is known and which has been used as a benchmark problem 
by others. 
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The problem involves the steady flow of a perfect gas in a rectangular 
duct in which density, velocity, and energy are prescribed in each of four 
triangular wedges in such a way that the appropriate jump conditions (the 
Rankine-Hugoniot conditions) are exactly satisfied. Thus, a problem of 
shock reflection for which an exact solution is known is obtained. Dimen- 
sions and data are given in Fig. 2. In this and all the other problems, 
the solution is considered to have converged to steady state when the 
magnitude of the L 2 -norm of the density is reduced by three orders of 
magnitude. 

The time step is monitored by the formula 



2 yP I 1 2 2 -'2 

Here, C = — and |u| = u^ + , Y = 1.40 . The constants 

multiplying the artificial viscous terms were selected locally as: 



where the bar denotes average element values. A Lapidus constant of 
1.0 is used. 

The results of a uniform coarse grid approximation are shown in 
Fig. 3. The computed density contours are also shown. 

The same problem was also analyzed using the node redistribution 
algorithm with 20 node redistribution iteration. Results are shown 
in Fig. 4. There, the original coarse initial mesh of Fig. 3 is progressively 
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Figure 2 . 


A shock reflection problem, 
variables are prescribed as 
outflow values are computed 
laws . 


Inflow values of the conservation 
indicated in regions I and II, and 
in III to satisfy the conservation 
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distorted to conform to the reflected shock locations. Corresponding 
density contours are also given in the figure. 

NACA 0012 Airfoil in Supersonic Wind Tunnel . In this example, 
the supersonic flow through a narrow wind tunnel containing a NACA 0012 
airfoil is studied. The inflow Mach number was set at M^ = 2 , with 
Y => 1.40 and symmetry is exploited to reduce the computational effort. 

The initial coarse mesh and density computed contours are given in 
Fig. 5. Note that the critical features of the solution — the 
reflected shock and contact discontinuity — are lost with this coarse 
mesh. Results of a node-redistribution scheme for the coarse mesh 
are shown in Fig. 6. In these results, ten iterations of the node re- 
distribution algorithm were used. 

Supersonic Flow in a Wind Tunnel with a Step . The steady-state 
solution of the problem of a wind tunnel with a step introduced Into the 
flow is next considered. The inflow Mach number was selected M^ = 3.0 
and Y = 1.40 . The initial coarse mesh is shown in Figure 14 with the 
corresponding density profiles. The mesh refinement algorithm was also 
used, with the mesh and density profiles obtained after 10 iterations 
shown in Fig. 8. We see that some oscillations are present downstream, 
and they are believed to be due to the non-monotonicity of the solution 
algorithm. The results presented for the refinement-unrefinement procedure 
have been constrained by a maximum number of 2000 nodes or 2000 elements 
that can be allowed. In the refined mesh shown, this constraint has 
been achieved . 
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Figure 3. Reflecting shock problem 


Initial mesh and density contours. 
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Figure 4 . Reflecting shock problem. Mesh and density contours obtained 
after 2x10 applications of the mesh redistribution algorithm. 
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Figure 5 . NACA 0012 airfoil in supersonic wind tunnel. 
Initial mesh and density contours. 
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Mesh and density 
applications of 


in supersonic wind tunnel. 
Dntours obtained after 2x10 


mesh redistribution algorithm. 




Figure 7 * Supersonic flow in a wind tunnel with a 

step. Initial mesh and density contours. 
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Figure 8 . Supersonic flow in a wind tunnel with a step. 

Mesh and density contours obtained after 2x10 
applications of the mesh redistribution algorithm. 
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9 . Blunt 1< 
Initial 



ag edge in hypersonic Elow field 
h and density contours. 
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Figure 10 . Blunt leading edge in 
hypersonic flow. Mesh 
and density contours 
obtained after 4 
applications of 
the mesh 
redistribution 





algorithm. 



Blunt Leading Edge of 8* HTT Panel Holder in Hypersonic Flow . 
The problem of the blunt leading edge of the 8' HTT panel holder In a 
supersonic flow field with freestream Mach number = 6.57 , y = 1.38 
and 0° angle of attack, was solved to obtain the steady-state solution. 

This problem has also been studied by Bey et al . 

A coarse mesh solution is indicated in Fig. 9. A distorted mesh 
and corresponding density map are indicated in Fig. 10. In this particular 
problem, the r-method did not give particularly good results, as a poor 
approximation of the solution between the shock and blunt body results 
from spurious ocsillations in the basic time-marching algorithm. In 
the case of mesh adaptation using redistribution, the solution actually 
diverges after four passes through the adaptive scheme due to the badly 
graded (hourglassed) mesh produced from the oscillations of the adaptive 
scheme downstream of the shock. 

3.4. A p-Method . Returning now to the full, viscous, Navier-Stokes 
problem characterized by (2.1), we outline a p-method for adaptive 
improvement of finite element solutions. The procedure is straightforward: 

1 . On a fixed mesh, compute a trial solution at each time step 

using a linear or bilinear approximation of the velocity field U^. 

£ 

2. Use to compute a residual for data in a local problem which 

yields the local error indicator <J> over each element K. 

K 

3. Use local quadratic approximations to compute approximate error 

indicators <p . 

K 
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4. Check for pre-assigned error tolerances. If the error is too 
large, recompute locally using quadratic elements. 

We employ hierarchal families of elements and could go up to quadic 
polynomials. In results given below, only quadratics are employed, 
and mesh refinement is used if increasing p fails to bring error with 
acceptable limits. 

To Lest the residual estimation technique and our p-method, we 
examine the performance of the method using a problem with a known 
analytical solution. Since few analytical solutions to the transient 
Navier-Stokes equations are known, we attempt to construct solutions 
to special problems by choosing an analytic solution and then computing 
the corresponding boundary conditions and body force. 

The theoretical analysis suggests the following local error indicators: 


{pU K |K,t + P 


. K , 

'dt - 1 ,K .t 


! }dt + 
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Here , 


ifl^t ( ! K *! K) t ; l^lx.t = (( + K » ! K)) t 


For each error indicator e,. , j 

k. , n 

effectivity ratio 


= 1,...,4 we define the 


0 j = _fLljL 

K,n 


' ?h'' K, n 


where lie II ^ is the norm estimated by ep, as described above, 
-h K, n J K,n 

say that is an effective estimator if the corresponding 

K, n 

effectivity ratio takes values close to unity. 


We 


In the numerical example to follow space-time elements of 

first order (bilinear in space, linear in time) are used to obtain an 

approximate solution and second-order space-time elements (biquadratic 

in space, quadratic in time) are employed to solve the local parabolic 

problems in the first time step. Then a normalized value ^ of 

K ,n 

the local error indicator (ffi is computed for every element K . A 
new order of approximation for element K is defined as follows: 


If 

0 « 3.. < 

6 

first order approximation 

If 

6 < V. ‘ 

1 

second order approximation 


Thus, the order of approximation is increased in the elements with 
error indicator bigger or equal to 6 times the value of the largest 
error indicator and an enriched mesh is obtained. This enriched mesh 
is used to advance the solution for M time steps (including the 
first) and after that the whole procedure is repeated. The success of 
this strategy in any particular problem depends on the choice of the 
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discretization parameters h , At and the constants 6 » M . 

We now consider a problem, with a one-dimensional "wave-like" 
solution in the fixed spatial domain Q = [0,20] x [0,20], n a m ely 

u = 0, v = 10e - ^ 0 * 8 ( x -4)“ 2t}2 , p = 0 

The problem was solved with two uniform (5x5 and 10 x 10) meshes 
using the following parameters: 

Time step: At = 0.025 , 0.100 

Penalty parameter: c = 10 3 

Fluid constants: p =1, y a 1 

Kinematic boundary conditions were applied on the side with y = 0 
while traction boundary conditions were applied on the remaining three 
sides of the boundary. This explains why the numerical solution" for 
the v-component shown in Fig. 11, 12 is not uniform in the y-direction. 

For each error indicator we may compute a global effectivity 
ratio defined by: 



Table 1 gives the evolution of the global effectivity ratios for 
the proposed error estimators through the first nine time steps. It 
is observed that the values of these ratios diverge from the optimal 
value of one as the time marching progresses. This phenomenon is due 
to accumulation of the truncation error which destroys the effectivity 
of the error estimates and makes the adaptive procedure inefficient as 
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X 


Figure 11. Problem with exponential solution. Computed solution for 
time t = 0.1. Parameters: h = 2, A t = 0.1, 6 = 0.25, 

M = 4 . 
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Fi gure 12.. 


Problem with 

time t = 0.9. 


e *ponential solution. 
Parameters: h = 2 


Computed solution f„ r 
° °‘ J . 6 - 0.25, Ms 


4. 
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can be seen. This situation can be corrected by reducing the values 
of the parameters of the adaptive procedure, namely h. At, M, 6. 


Table 3.1 

Values of the global effectivity ratios after the first time step; At 
0.025 


Mesh 

e i 



ef 

5x5 

1.0439 

1.6995 

1.3037 

1.0339 

10x10 

1 . 0054 

0.3225 

1.1627 

1.0011 
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4. CONCLUSIONS 


1. The use of rigorous residual error estimates for time-dependent 
Navier-Stokes problems in two dimensions is justified by the methods and 
results developed during the course of the work. In particular, sharp 
error bounds can be derived which were found to yield excellent error 
estimates for a test problem. 

2. The use of p-methods and residual error estimates can be computa- 
tionally expensive. New data structures must be developed before these 
types of methods can find wide application in practical flow problems. 

3. Adaptive p-methods may ultimately offer significant advantages 
over other adaptive schemes if used in parallel computation with efficient 
use of multiple array processes. At present, these schemes do not perform 
well in comparison with some other adaptive strategies from the viewpoint 
of programming ease and computational efficiency. 

4. Node redistribution methods are simple to use and can be effective 
in certain steady problems. They do not appear to be attractive in general 
flow problems, particularly in cases in which boundary conditions vary in 
time . 

5. The node-redistribution-relaxation scheme described in this report 
can exhibit numerical instabilities in certain flow problems. The method 
is incapable of reducing local errors below very high tolerance levels 
without acquiring some destabilizing properties. It is doubtful that the 
method will find broad applications unless it is combined with an h-type 
or p-type adaptive strategy. 
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6. The h-methods, not dealt with here, do offer some advantages in 
terms of speed of implementation and general data structures over the 
r-method and the p-method developed here. However, these schemes can also 
exhibit instabilities. It is highly likely that the best adaptive schemes 
will prove to be those which combine h- and r-strategies . 

7. Interpolation methods for error estimates are rather easy to 
implement and can be easily added to existing CFD software. They can 
provide only a crude indication of the actual local errors, but they can 
be used to produce an adequate indication of relative error between succes- 
sive meshes. 

8. Additional studies are needed to explore algorithms for accurate 
and efficient time-dependent problems. While we have developed some schemes 
in earlier work (e.g. [1]), a general approach to this important problem 
which fulfills requirements of accuracy and efficiency is not yet available. 
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APPENDIX 


A Posteriori Error Estimation for a Class of One-Dimensional Elliptic 
Problems. 

Consider the model problem, 
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( 1 ) 


( a ■“) + bu = f in (0,1) ] 

dx dx ! 


u(0) = u(l) = 0 


where it is assumed that 


0 < a^ < a(x) < a^ < + ~ Vxe [0,1] 


0^b(x)<b 1 Vxe [0,1]. 

Assuming a,a',b,feC([0,l]) we have ueC 2 ([0,l]). 


Let us define the bilinear form, 


- 1 du dv 


a(u,v) = f (a + buv)dx 

J dx dx 


Then u is the solution of the variational equation, 
l 

a(u,v) = j fvdx VveHj ([0,1]) 

0 

Moreover, we define the energy norm ||-||g by. 


IMIe = a(v,v) VveH^O.l]) 

We discretize [0,1] using N linear element Qj = [xj_j, xj], j=l,...,N, where 
0=x 0 <X[< • • • <Xj_]<Xj • - ' <X]sj=l. Let V h cHo ([0,1]) denote the space which is 
spanned by the basis functions of the discretization. We approximate the solu- 
tion u of (1) with a u h €V h which is obtained as the solution of the discrete 
problem: 
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Find u h e V h such that, 


a K> v h) = J fv hd* Vv h e yh 

o 

We want to estimate the energy norm of the approximation error, 
e h = u-u h . In view of the orthogonality condition, 

a(e h , v h) = ° Vv h eV h 

the energy norm of the approximation error satisfies, 

IMe = a ( e h- e h - v h) 

* u * A 

Choose v h eV" such that u h (Xj)=e h (xj), j=0,l,...,N and let w=e h -u h to obtain: 

IMe = E T ( a ^(u-u h )-^- + b(u-u h )w}dx (2 ) 

j=0 ^ dx dx 

Integrating by parts and noting that u^=0 in (xj, Xj+j), j=0,l,...N-l, we get, 

l|e h lli = X T (f+aV h -bu h )wdx (3) 

j=0 Xj 

We also have (see Babuska and Rheinboldt [i]) that: 

H e h- v hllE * (l+0(h))||e h || E (4) 

Starting from (2) and using Schwarz’s and Korn’s inequality and (3) we get 

the following estimate: 
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( 5 ) 


||e h || E <(l+0(h))K(h)( N f V / (f+a'u' h -bu h ) 2 dx)* 

J=0 Xj 

Here, h = max (x:. i-X:), and K(h) denotes Korn’s constant. 

0<i<N-l 3 3 

The estimate given in (4) has the disadvantage that it involves Korn’s 
constant which has to be computed for a given discretization and even then the 
estimate is far from being optimal. Babuska and Rheinboldt [2], developed a 
less obvious approach ( than the one leading to (4) ) which gives sharp error 
estimates. This approach is summarized below. 

Define a family of local auxiliary problems as follows: 

For every element Qj +1 , find 4)j +! e H 0 1 ([x J , Xj +1 ]) such that, 

a (<t>j+i. v ) = J (f+aV h -bu h )vdx 
VveHjdxj, x j+1 ]), j=0,...,N-l 

Combining (2./- with (5) using Schwarz’s inequality and inequality (3 ) we 
obtain the bound: 

Me * (2 

i=0 

One can estimate ||<t>j + illE,n H using the Lax-Milgram theorem to obtain: 

H4 > j+i!lE,n H - ^j+iH f+a,u, h -bu hllo,fVi 

Moreover by using spectral analysis one can show (see [ 3 ]) 
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( 6 ) 


(x j+ 1 -Xj ) 2 

H+i s 2 ■ TT 

min a(x) 

xe(Xj^„ 

and finally we get the estimate: 


' 'h 

(7) 

We note that the error bound given in (6) can be computed without solving 
the local problems defined in (5). 


INIe ^ 


N £ ( x i+i~ x j) 

j^o 7t 2 min a(x) 

X£(Xj^]) 


||f+a'u' h -bu h || 0 2 ^ : 


•H 


Locally Computed A-Posteriori Error Estimates for Elliptic Equations 

For the sake of simplicity, we first restrict ourselves to the following 
model elliptic problem. 

Find u such that 


Au = f in Q 


u = u 0 on T u 

9u _ 

— = on It 
on 

Here Q is a bounded domain in R 2 , T u , and T t form two disjoint parts of the 
boundary dQ and f, u 0 , g are given functions. 

If u 0 denotes a function from H*(Q) whose extension to I~ u coincides 
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with u 0 ( when we use the same symbol for u<) and its extension) we have the 
usual variational formulation to ( 7 ) . 

Find u e uq + V such that 

JVu-Vv d£2 = j f-v dn + J gv ds V v e V 
n n r T 

where 

V = {u e = 0 on T u ) 

and d Q, ds denote measures of volume and surface area’respectively. 

In the following, we discuss a method to compute a-posteriori error esti- 
mates for finite element solutions of (7) . 

Suppose we are given two approximation spaces, say V h l , c V such 
that approximation (interpolation particularly) of an arbitrary element v e V by 
elements vf e Vf results in "much smaller" error than the approximation by 
vq 1 e V h \ Specifically we can think about V|P as the space spanned by the 
hierarchical elements with order p large enough. V,} does not mean necessarily 
the space corresponding to the first order approximation. 

The approximate problem for the "lower order" approximation can be 
stated as follows: 

Find ujj e uq + V h ] such that 
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jYuh-Vvh 1 dQ = J fvj dtt + jgVh 1 ds v Vh 1 € 
n n r T 

For simplicity, we will assume that kinematic boundary conditions can be 
satisfied exactly by the first order approximation (i.e., u 0 e Vj}). 

The residual rj corresponding to the solution u h ! is an element from the 
dual space V* and should be estimated in the corresponding dual space norm: 

ll r h llv = sup <r h 1 , v> (8) 

V 6 V 
IMIv*! 

Certainly r h is weakly continuous over V and the ball ||v|| v < 1 is weakly com- 
pact, and, therefore, there must exist a certain v 0 , ||v 0 |lv ^ 1 in fact it must be 
normalized in the sense that ||v 0 || =1) such that 

IMv = v o> 

* 

Now, let vR denote the best approximation of v 0 out of the space Vf. Suppose 
||v 0 - v£|| < e. Defining v£ = v£ /||v^|| we easily conclude that also 

II Vo - v&H ^ e + 

which is equivalent to saying that v 0 can be approximated "close enough” by ele- 
ments with unit norm. 

Next, ( 8 ) can be replaced by 
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( 9 ). 


Ilihllv* £ ll r hllv*H v (r v ^ll + sup <r£, v^> 

llvfll s 1 

the first term on the right-hand side being negligible. In other words we 
presume that the error in terms of the dual space norm can be estimated 
sufficiently well by taking the supremum over only the subspace Vjf. 

Now let us try to estimate this supremum. We have: 

<Jh> v £> = jYuh'Vv,? dn-JfvgdQ- ( T gvf ds = 

n n 

, . i 9uh 5uu* 

- Au h J - f)vg dfl + j -(— T— K ds + 

aK-3£2 2 ^ 

r ^ u h 

+ f (— g)V(? ds (2.n+y) 

3K n r T dn 

Let us note also that in (2.24) v£ can be replaced by the diffemce v£ - vj}, with 
arbitrary v,} e VjJ since the residual r h l is orthogonal to V,j. In particular, we 
can choose for - a Vj - inteipolant of v|. 

Finally, defining 

V&(K) = (v£ e Vp)|V‘ - interpolant of vjf = 0} 
we may use the following local problem associated with an element K: 

Find 4> K € VjP(K) such that: 


K K 
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jY4>icY v f? = 

K 


|(-Vu,}— f)V(f d£2 + / 

k 3K-an 


1 (ilii 

2 C dn 


3uh U 

dn 


)V[P ds 


r 

+ | (-T- - g)vjP ds V vfo € VjP(K) 

inr T 30 

and equation ( 9 ) can be rewritten in the form: 


<Th> v £> = £ JY<1 >kY(V|? - v,, 1 ) dQ 
K k 

The functions <J>k are local error indicator functions. 

Finally, we note that for each element K there exists a constant C K 
such that 

W - Vj'lb.K S C K ||vf||,, K V yg € Vf(K) 
which results in 

«h‘. vff> S C(£ J(Y« 2 dn^llvjllv 

K K 

provided 


C = maxCif 

K K 

The constant is independent of the size of element K, but it does depend on 
its shape. In particular, for grids consisting of the same regular elements (not 
necessarily uniform C K is independent of K and does not affect the local 
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behavior of the estimate (10) (comp, also [ 4 ]). 


Our estimate: 

llr h V * C(£ J(Y * k ) 2 dH)* 

KK 

clearly has only a global character; however, the contributions J(Y4 >k) 2 dO 

K 

(local) error indicators) may serve as a basis for some local refinement 

Finally, let us note that in the case of the model problem the residual 
estimate can be directly reinterpreted as the estimate of the error ej = u h ’ - u, 
since 


!l e h llv = |Y e h Y e h dH = 


n 


= JVuhVeJ d£i - J(-Au)eh dO - J ds = 


n n 

<Th> eh> ^ 11^11*11^11 
and therefore. 


5n 


lien 1 llv ^ ll r hll* 

Conversely taking into account that in the case of a Hilbert space H with inner 
product (.,.) H and the corresponding norm ||-|| H , 


||u|| H = su 
ve 


<U, V> 

8 IMI 


one gets that 
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ih‘n. = iK’iiv do 

In fact, for any self-adjoint elliptic problem we may define the proper 
(energy) norm in such a way that the norm of residual is exactly equal to the 
norm of the error. For a class of not self-adjoint problems or for certain classes 
of nonlinear operators (strongly monotone) the residual norm can likewise be 
bounded below by the norm of the error. However, in general the relation 
between these two quantities is more complicated and we must be satisfied with 
the observation that, in general when the solution of the boundary value problem 
depends continuously on the righ-hand side (the data), it makes sense to minim- 
ize the residual since it results in simultaneous convergence of the approximate 
solution to the exact one (see also [ 5 }). 

We now show that the concept of estimating the residual over a higher 
order approximation space can be easily applied to Stokes problem. Stokes 
problem is of special interest to us since it is obtained from a linearization of the 
stationary Navier-Stokes equations. 

Let us restrict ourselves to the case of pure kinematic boundary condi- 
tions (the Dirichlet problem). The Stokes consists in the determination of the 
velocity field u and pressure field p such that: 
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- jxAu + grad p = f 

div u = 0 in Q 
u = u 0 on dQ 

where (I is the viscosity constant, f stands for body forces and u 0 specifies velo- 
city on the boundary. The usual consistency condtion must be satisfied, i.e., 

J u,j ds = 0 •* 
an 

where u n = un, provided n denotes the outward normal unit to dQ. The varia- 
tional formulation of (2.32) is as follows: 

Find u 6 u 0 + V, p e Lq(D) such that: 

\ij e ij (u)e ij (y) dD + J p div v dQ = J fv dD V v e y 
n n n 

Jdiv u q d£i = 0 V q €^Lq(Q) 
n 

Here e,j(u) denotes the symmetric part of the gradient of u. In the above u 0 
denotes an arbitrary function in such that uqI^q = uq (again, we use the 

same notation for u 0 and for its trace on the boundary, and that div u 0 = 0, and 
Y = Hq(^) and, as usual, 

Lo(Q) = {q e L ? (D)| J q dQ = 0} 
n 

Denoting, as previously, a p-th order approximation space by Yf and 
Lo h p respectively, we may formulate the p-th order approximate problem as fol- 
lows: 
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Find ufi € + Vfi, pfi e L^ h (ft), such that 


pj ^(ufJeijCv^) dQ. + /p^P div (yJP) dT2= J f yg dQ. V vjf € V;P ( 12 ) 
n n n 

J div q£ dfl = 0 V q e 4 1 ! 
n 

As in the preceding case, we attempt to estimate the proper residuals 
corresponding to the "1-st order approximation’ , i.e., replacing suprema over 
v € V and q e Lq(Q) by suprema over vfi e Vfi and qj, e L^jf, with some 

p > 1. 

We have for the first equation: 

< “ + grad ph 1 - f; yf > = 

Z^jCUhteijCyi?) dfl ~ JPh div v| dn - j fv£ dD} = 

K K K K 

Z(J(- Aa h + grad Ph 1 - f)v£ d£2 + 

K k 

J • t(y>1 ' 1> “ l) -- 1 * ( “ l ’ 1,P “) . y| P ds } (13) 

3K an 2 

where, as usual, jt* denotes a stress vector corresponding to an adjacent element. 
For the second equation we get simply: 

< div u,, 1 , qj? > = £ J div u,, 1 qj? dn 

Kk 

As usual, one can see that vf and q£ in ( 12 > and ( 13 ) can be replaced by 
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v£ - v h ' and qjf - qjJ for arbitrary V), 1 and particularly choosing for v^ and 
interpolants of y£ and q£ we can set the local problem in the following way: 

Find $ K e Y^ 0 (K), y K e L&(K) such * at: 

pj(- Aud grad ph 1 - f) yf P dQ + j ^ygdQ. Vy^ £ Y£o(K) 

K 8K XI 1 

j div($ K )q£ dQ = 0 V q£ e L&CK) 

K 

In the above Y£,(K) = {u| e YifCK^Yh'^^T 50 ! 2111 = 0) 

L£(K) = {q£ € L^K)^ 1 - interpolant of q^ = 0} 

Finally (.1-2 ) and (13) can be rewritten in the form: 

< - pAujJ + grad - f,vfi > = 

X (nj £ij(^K)e.j(y^-yh) dfl - JVk div (y^-Yh 1 ) d^ 

K K K 

and 

< div Uh 1 , q^ > = Xj divC^^CqjP-q^ 5 dQ 

Kk 

where vj and qjJ are the appropriate interpolants of and q£ respectively. Both 
identities lead to the residual error estimates. However, the second one is practi- 
cally useless, since the ||div ujJ ll L *( 0 ) can be estimated simply by a direct compu- 
tation. Let us therefore focus our attention on the first residual. By Schwarz’s 
inequality, we get 
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Inj e 1J (^>K) £ 1 j(y^-yh 1 ) dOMEGA - Jy K div (y^-y,, 1 ) dfi| £ 

K K 

^ (II4>kIIe,k + II VkIIojc)^II V J — v h IIk ^ 

^ CkOI^He.k + HVkIIo.kJ^II^Hk- 

where ||*|Ih,k is energy norm associated with the viscous, bilinear form, 
||v|| K = (Ilyll^K + lldiv y||o >K ) v ’ is a well-defined norm in V^ 0 (K), and finally C K is 
a constant associated with the estimate of the type: 

IK - yiillK s CkIWIIk v yjf £ v,f(K) 

C K is independent of the size of the element K but depends upon its shape. For 
detailed discussion of these ideas, we refer to [M • We note here that for 
meshes which are uniform or consist of elements of the same shape the constant 
C K is independent of the element and does not affect the local behavior of the 
estimate. 

Finally, we arrive at the estimate: 

II - + gradPh ~/ll-i,n - CII^kIIe.k + IIVkIIo,k)) /j 

K 

where C = max C K , and the norm in the dual space is measured with respect to 
K 

the norm: 


||v|| = (I (1 |v|£ k + lid* v IIo,k))^ 

K 
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